

%center
x2 = 0; 
y2 = 0;

%location of receiver A
x1 = 0; 
y1 = 10; 

%location of receiver C
x3 = 10; 
y3 = 0;

%epicenter
epicenter = [7 3];


R1 = sqrt((epicenter(1) - x1)^2+(epicenter(2) - y1)^2);
R2 = sqrt((epicenter(1) - x2)^2+(epicenter(2) - y2)^2); 
R3 = sqrt((epicenter(1) - x3)^2+(epicenter(2) - y3)^2); 

%recordered TDOA
R12 = R1 - R2 ;
R21 = -R12;
R32 = R3 - R2;
R23 = -R32;
R13 =  R1 - R3;
R31 = -R13;

X = 5-1/20/(4*R32^2+4*R12^2-400)*(400*R12-4*R12^3+8*R32*R12^2+4*R32^3-400*R32-4*(2000000-200*R12^3*R32+400*R32^2*R12^2+100*R12^4+100*R32^4-R12^4*R32^2+2*R12^3*R32^3-200*R12*R32^3-R32^4*R12^2+20000*R12*R32-30000*R32^2-30000*R12^2)^(1/2))*R32+1/20*R32^2;

%X = 1/2/(400*R12^2+400*R32^2-40000)*(-40*R32^2*R12^2+4000*R12^2-400000+4000*R32^2-4000*R12*R32+40*R12^3*R32+40*(100*R32^6-30000*R32^2*R12^2-30000*R32^4+100*R12^4*R32^2-200*R12^3*R32^3+20000*R12*R32^3-R32^4*R12^4+400*R32^4*R12^2+2000000*R32^2-200*R12*R32^5+2*R12^3*R32^5-R32^6*R12^2)^(1/2)); 

Y = 5-1/20/(4*R32^2+4*R12^2-400)*(400*R12-4*R12^3+8*R32*R12^2+4*R32^3-400*R32-4*(2000000-200*R12^3*R32+400*R32^2*R12^2+100*R12^4+100*R32^4-R12^4*R32^2+2*R12^3*R32^3-200*R12*R32^3-R32^4*R12^2+20000*R12*R32-30000*R32^2-30000*R12^2)^(1/2))*R12-1/20*R12^2+1/10*R12*R32;
 
% Y = 1/20*(100*R32+10*R12/(400*R12^2+400*R32^2-40000)*(-40*R32^2*R12^2+4000*R12^2-400000+4000*R32^2-4000*R12*R32+40*R12^3*R32+40*(100*R32^6-30000*R32^2*R12^2-30000*R32^4+100*R12^4*R32^2-200*R12^3*R32^3+20000*R12*R32^3-R32^4*R12^4+400*R32^4*R12^2+2000000*R32^2-200*R12*R32^5+2*R12^3*R32^5-R32^6*R12^2)^(1/2))-100*R12-R32*R12^2+R32^2*R12)/R32;

figure(4)
plot(X, Y, 'x'); 
